// Program to compute non-parametric LLR (ridge) regression in MATA
mata:
real matrix llr(real colvector p0, real colvector  p1, real colvector  y0, real rowvector h, real rowvector kernel)
{
   real matrix dist, G
   real scalar zero, n, c
   real colvector s0, miss, x1, distx1, s2, t0, t1, ridge, yout
   dist = (p1 :- J(rows(p1), cols(p0),1) * p0')'
   zero = 2.22045e-16   //See MATA Manual, Page 320
   n=rows(p0)
   if (h==10000) dist = J(rows(p0), rows(p1), 0)
   if (kernel==1) G = ( abs(dist) :< h ) :* ( 1 :- (dist:^2):/ (h:^2 ))*(3/4) //epanechnikov
   if (kernel==2) G =  normalden(dist:/h)                                     //gaussian
   if (kernel==1) c=0.3125      //epanechnikov
   if (kernel==2) c=0.3535      //gaussian
   //s0  = colsum(G:* dist:^0 )'
   s0  = colsum(G)'
   miss = abs(s0):<=zero
   x1  = editmissing((1:-miss):*(colsum(G :* p0)' :/ s0),0) :+ miss:*p1
   distx1= (p0 :- J(rows(p0), cols(x1),1) * x1')
   s2  = colsum(G:* distx1:^2 )' // s1 is zero by construction of x1
   t0  = colsum(G:*y0)'
   t1  = colsum(G:*distx1:*y0)'
   if (h>0 & h<10000) ridge=h:*abs(p1-x1):*c
   if (h==10000)      ridge=999999
   yout = editmissing((1:-miss):*(t0:/s0 + (p1-x1):*t1:/(s2:+ridge)),0)
   yout = editvalue(yout, 0, .)
return(yout)
}
end

// Program to do cross-validation
mata:
real matrix cv(real colvector p0, real colvector p1, real colvector y0, real colvector y1, real rowvector kernel)
{
   real matrix R, vec, hopt
   real scalar i,j, g, h, n0, N, mincv, minbound
   real colvector hvalues, cv, yhat0, pp0, yy0,  yhat, mark, y
   //Grid
   hvalues = J(30,1,0);
   hvalues[30,1] = 10000;
   for (g=1; g<=29; g++){
     hvalues[g,1] = (0.01)*(1.2)^(g-1)
   }
   // Initialize cv function   
   cv = J(rows(hvalues), 1, .)
   // Compute yhat
   for (j=1; j<=rows(hvalues); j++){
    h = hvalues[j,1]
    n0    = rows(p0)
    yhat0 = J(n0,1,.)
    pp0   = p0[(2..n0),.]
    yy0   = y0[(2..n0),.]
    for (i=1; i<=n0-1; i++){
     yhat0[i,1] = editmissing(llr(pp0, p0[i,1], yy0, h, kernel),0)
     pp0[i,1]  = p0[i,1]
     yy0[i,1]  = y0[i,1]
    }
    yhat0[n0,1] = editmissing(llr(pp0, p0[n0,1], yy0, h, kernel),0)
    y= y0
    yhat = yhat0
    N = n0
    cv[j,1]   = (1/N)*colsum( (y-yhat):^2)
   }
   //Find minimum
   cv = editvalue(cv,0,.)
   mincv     = min(cv)               // internal min //
   minbound  = cv[rows(hvalues),1]   // min on the boundary //
   if (minbound>mincv) {; mark  = cv:==mincv;   };
   if (minbound<=mincv){; mark  = cv:==minbound;};
   vec   = hvalues:*mark
   hopt  = J(rows(hvalues),1,1):*max(vec)
   R = (hopt,hvalues, cv)
return(R)
}
end

// Program to compute LLR (ridge) treatment effect estimator in STATA
capture program drop llrmatch
program define llrmatch, rclass
qui {
syntax varlist(max=2) [if], Pscore(string) Kernel(string)
** Parse Data
   tokenize `varlist'
   local Y = "`1'"
   local T = "`2'"
   local pX="`pscore'"
   tempvar d0 d1 sort_preserve sample
   gen `sample' = 0
   replace `sample' = 1 `if'
   gen double `d0' = `T'==0 & `sample'==1
   gen double `d1' = `T'==1 & `sample'==1
   mata: p0 = st_data(.,"`pX'","`d0'")
   mata: p1 = st_data(.,"`pX'","`d1'")
   mata: y0 = st_data(.,"`Y'","`d0'")
   mata: y1 = st_data(.,"`Y'","`d1'")
** Control kernel type option
   if "`kernel'"=="epan"     local kN = 1
   if "`kernel'"=="normal"   local kN = 2
** Cross-validation
   mata: hopt = cv(p0,p1,y0,y1,`kN')
   mata: hstar = mean(hopt[.,1])
   matrix define hstar = . 
   mata: st_replacematrix("hstar", hstar)
   local hstar = el(hstar,1,1)
** Match and gen weight variable if outcomes not specified;
  mata: Y0hat = llr(p0, p1, y0, `hstar', `kN')
  mata: tot=mean(y1:-Y0hat)
  mata: st_numscalar("tot",tot)
  return scalar llr = tot
}
end


